R jako GIS

Podpůrný text pro předmět: GIS pro biologické aplikace
Autor: Matěj Man
Aktualizace: 08. 11. 2021
Odkaz na R skript: Kopírovat a vložit do R
Odkaz na data: Stáhnout data zip Odkaz na github repozirář: Repozitář

Další zdroje:

Knihovny

# list.of.packages <- c("sf","raster","mapview","whitebox","randomcoloR","leaflet","Rcpp")
# install.packages(list.of.packages)
# install.packages("devtools")
# install.packages('raster', repos='https://rspatial.r-universe.dev')
# library(devtools)
# install_github("r-spatial/sf")

library(sf)
library(raster)
library(mapview)
library(randomcoloR)
library(leaflet)

Načítání vektorových GIS dat do R

## Načítání vektorových dat shp
# nastavte kde máte u sebe na PC data
# pozor musíte zdvojit nebo otočit lomítka
cesta<-"d:/Git/gisproba/data/"
# zkonstuuje cestu kde leží vrstva Brdy
data.path<-paste0(cesta,"CHKO_Brdy.shp")
# načte .shp do R
brdy<-st_read(data.path,stringsAsFactors = F,quiet = T)

# obdobně načteme třeba hranici ČR
data.path<-paste0(cesta,"hrcr1_wgs.shp")
hrcr<-st_read(data.path,stringsAsFactors = F,quiet = T)

# prostý obrázek, bez interaktivity
plot(st_geometry(hrcr)) 
# parametr add přidá další vrstvu do existujícího obrázku
plot(st_geometry(brdy),col="red",add=T) 

# Nastaví pracovní adresář working dorectory (WD)
# dále nemusíme data z WD volat celou cestou, stačí názvy 
setwd(cesta) 
## Načítání vektorových dat z tabulky
# načíst prostou tabulku
chmu<-read.table("staniceCHMUtablecoma.csv",header = T, sep=",")
# převést tabulky na prostorová data
chmu_sf<-st_as_sf(chmu, coords = c("Xcoo", "Ycoo"),crs = 4326) 

## vizualizovat prostorová data interaktivně
# pokročilejší balíček leaflet

# definujeme paletu o 17 barvách podle typu stanice
distCol<- colorFactor(distinctColorPalette(17), chmu_sf$Station.ty)

# definoujeme okno s mapou 
leaflet(chmu_sf) %>% 
  addTiles() %>%  
  addCircleMarkers(popup=chmu_sf$Name,color = ~distCol(Station.ty),
                   stroke = FALSE, fillOpacity = 0.8,radius=4) %>%
  addLegend(pal = distCol, values = ~Station.ty)

Projekce a crs

st_crs(brdy) # jaký crs má vrstva brdy ?
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
st_crs(hrcr) # jaký crs má vrstva hrnice čr ?
## Coordinate Reference System:
##   User input: WGS 84 + Unknown VCS from ArcInfo Workstation 
##   wkt:
## COMPOUNDCRS["WGS 84 + Unknown VCS from ArcInfo Workstation",
##     GEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]],
##         CS[ellipsoidal,2],
##             AXIS["longitude",east,
##                 ORDER[1],
##                 ANGLEUNIT["Degree",0.0174532925199433]],
##             AXIS["latitude",north,
##                 ORDER[2],
##                 ANGLEUNIT["Degree",0.0174532925199433]]],
##     VERTCRS["Unknown VCS from ArcInfo Workstation",
##         VDATUM["Unknown"],
##         CS[vertical,1],
##             AXIS["gravity-related height (H)",up,
##                 LENGTHUNIT["metre",1,
##                     ID["EPSG",9001]]]]]
## vektory
# transformace podle EPSG
brdy_32633<-st_transform(brdy,32633)
st_crs(brdy_32633) # jaký crs má vrstva brdy_32633 ?
## Coordinate Reference System:
##   User input: EPSG:32633 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 33N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",15,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State."],
##         BBOX[0,12,84,18]],
##     ID["EPSG",32633]]
# transformace podle proj4 string
brdy_5514<-st_transform(brdy, "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")
st_crs(brdy_5514) # jaký crs má vrstva brdy_5514 ?
## Coordinate Reference System:
##   User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 
##   wkt:
## BOUNDCRS[
##     SOURCECRS[
##         PROJCRS["unknown",
##             BASEGEOGCRS["unknown",
##                 DATUM["Unknown based on Bessel 1841 ellipsoid",
##                     ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
##                         LENGTHUNIT["metre",1,
##                             ID["EPSG",9001]]]],
##                 PRIMEM["Greenwich",0,
##                     ANGLEUNIT["degree",0.0174532925199433,
##                         ID["EPSG",9122]]]],
##             CONVERSION["unknown",
##                 METHOD["Krovak (North Orientated)",
##                     ID["EPSG",1041]],
##                 PARAMETER["Latitude of projection centre",49.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8811]],
##                 PARAMETER["Longitude of origin",24.8333333333333,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8833]],
##                 PARAMETER["Co-latitude of cone axis",30.2881397222222,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",1036]],
##                 PARAMETER["Latitude of pseudo standard parallel",78.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8818]],
##                 PARAMETER["Scale factor on pseudo standard parallel",0.9999,
##                     SCALEUNIT["unity",1],
##                     ID["EPSG",8819]],
##                 PARAMETER["False easting",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8806]],
##                 PARAMETER["False northing",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8807]]],
##             CS[Cartesian,2],
##                 AXIS["(E)",east,
##                     ORDER[1],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]],
##                 AXIS["(N)",north,
##                     ORDER[2],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]]]],
##     TARGETCRS[
##         GEOGCRS["WGS 84",
##             DATUM["World Geodetic System 1984",
##                 ELLIPSOID["WGS 84",6378137,298.257223563,
##                     LENGTHUNIT["metre",1]]],
##             PRIMEM["Greenwich",0,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             CS[ellipsoidal,2],
##                 AXIS["latitude",north,
##                     ORDER[1],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##                 AXIS["longitude",east,
##                     ORDER[2],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##             ID["EPSG",4326]]],
##     ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
##         METHOD["Position Vector transformation (geog2D domain)",
##             ID["EPSG",9606]],
##         PARAMETER["X-axis translation",570.8,
##             ID["EPSG",8605]],
##         PARAMETER["Y-axis translation",85.7,
##             ID["EPSG",8606]],
##         PARAMETER["Z-axis translation",462.8,
##             ID["EPSG",8607]],
##         PARAMETER["X-axis rotation",4.998,
##             ID["EPSG",8608]],
##         PARAMETER["Y-axis rotation",1.587,
##             ID["EPSG",8609]],
##         PARAMETER["Z-axis rotation",5.261,
##             ID["EPSG",8610]],
##         PARAMETER["Scale difference",1.00000356,
##             ID["EPSG",8611]]]]
# transformace podle jiné existující vrstvy
brdy_krovak<-st_transform(brdy, st_crs(brdy_5514))
st_crs(brdy_krovak) # jaký crs má vrstva brdy_5514 ?
## Coordinate Reference System:
##   User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 
##   wkt:
## BOUNDCRS[
##     SOURCECRS[
##         PROJCRS["unknown",
##             BASEGEOGCRS["unknown",
##                 DATUM["Unknown based on Bessel 1841 ellipsoid",
##                     ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
##                         LENGTHUNIT["metre",1,
##                             ID["EPSG",9001]]]],
##                 PRIMEM["Greenwich",0,
##                     ANGLEUNIT["degree",0.0174532925199433,
##                         ID["EPSG",9122]]]],
##             CONVERSION["unknown",
##                 METHOD["Krovak (North Orientated)",
##                     ID["EPSG",1041]],
##                 PARAMETER["Latitude of projection centre",49.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8811]],
##                 PARAMETER["Longitude of origin",24.8333333333333,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8833]],
##                 PARAMETER["Co-latitude of cone axis",30.2881397222222,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",1036]],
##                 PARAMETER["Latitude of pseudo standard parallel",78.5,
##                     ANGLEUNIT["degree",0.0174532925199433],
##                     ID["EPSG",8818]],
##                 PARAMETER["Scale factor on pseudo standard parallel",0.9999,
##                     SCALEUNIT["unity",1],
##                     ID["EPSG",8819]],
##                 PARAMETER["False easting",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8806]],
##                 PARAMETER["False northing",0,
##                     LENGTHUNIT["metre",1],
##                     ID["EPSG",8807]]],
##             CS[Cartesian,2],
##                 AXIS["(E)",east,
##                     ORDER[1],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]],
##                 AXIS["(N)",north,
##                     ORDER[2],
##                     LENGTHUNIT["metre",1,
##                         ID["EPSG",9001]]]]],
##     TARGETCRS[
##         GEOGCRS["WGS 84",
##             DATUM["World Geodetic System 1984",
##                 ELLIPSOID["WGS 84",6378137,298.257223563,
##                     LENGTHUNIT["metre",1]]],
##             PRIMEM["Greenwich",0,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             CS[ellipsoidal,2],
##                 AXIS["latitude",north,
##                     ORDER[1],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##                 AXIS["longitude",east,
##                     ORDER[2],
##                     ANGLEUNIT["degree",0.0174532925199433]],
##             ID["EPSG",4326]]],
##     ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
##         METHOD["Position Vector transformation (geog2D domain)",
##             ID["EPSG",9606]],
##         PARAMETER["X-axis translation",570.8,
##             ID["EPSG",8605]],
##         PARAMETER["Y-axis translation",85.7,
##             ID["EPSG",8606]],
##         PARAMETER["Z-axis translation",462.8,
##             ID["EPSG",8607]],
##         PARAMETER["X-axis rotation",4.998,
##             ID["EPSG",8608]],
##         PARAMETER["Y-axis rotation",1.587,
##             ID["EPSG",8609]],
##         PARAMETER["Z-axis rotation",5.261,
##             ID["EPSG",8610]],
##         PARAMETER["Scale difference",1.00000356,
##             ID["EPSG",8611]]]]
# kde jsou ty brdy? no nevidím je protože to má jiný CRS
plot(st_geometry(hrcr),main="Kde jsou ty Brdy?") 
plot(st_geometry(brdy_32633),add=T) 

# musíme tedy transformovat jedno nebo druhé
hrcr_32633<-st_transform(hrcr, 32633)
plot(st_geometry(hrcr_32633),main="No jo, už máme správný CRS")
plot(st_geometry(brdy_32633),add=T,col="red")

Načítání rastrových GIS dat do R

# načítání jdnoduchou funkcí "raster"
DEM<-raster(paste0(cesta,"DEM_Jested_100m.tif"))
# kontrolí obrázek
# data tu jsou
plot(DEM)

# ale leží na správném místě? 
# ukáže crs vrstvy
DEM@crs
## CRS arguments:
##  +proj=krovak +axis=swu +lat_0=49.5 +lon_0=24.8333333333333 +alpha=0
## +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs
# zkušeným okem, nebo podle errorů při další alanýze  poznáte, 
# že nemá dobře definovaný crs, je  to křovák, ale šptně definovaný.
# tak mu řekneme, jaká je korektní definice křováka
# http://freegis.fsv.cvut.cz/gwiki/S-JTSK

DEM@crs<-crs("+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")

# i rastr je možné vizualizovat interaktivně
# definujeme barvy
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(DEM),
                    na.color = "transparent")

leaflet(DEM) %>% 
   addTiles() %>%  
   addRasterImage(DEM,colors = pal, opacity = 0.8) %>%
   addLegend(pal = pal, values = values(DEM),
             title = "Elevation [m]")

Příklad vektorové analýzy

# plocha polygonů
st_area(brdy)
## 343898692 [m^2]
# délka linií (obvod polygonů)
st_length(brdy)
## 0 [m]
# buffer 5 km
brdy_5000<-st_buffer(brdy_32633,5000)
plot(st_geometry(brdy_5000), main="buffer 5 km")
plot(st_geometry(brdy_32633),add=T,col="red")

# centroidy
brdy_cen<-st_centroid(brdy_32633) 
## Warning in st_centroid.sf(brdy_32633): st_centroid assumes attributes are
## constant over geometries of x
plot(st_geometry(brdy_32633), main="centroid",graticule=T,axes=T)
plot(st_geometry(brdy_cen),add=T,col="blue",pch=19)

# průnik
chmu_brdy<-st_intersection(brdy,chmu_sf)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(st_geometry(brdy), main="intersection - chmu stanice v brdech",graticule=T,axes=T)
plot(st_geometry(chmu_brdy),add=T,col="red",pch=19)

Hromadné zpracování dat

Aneb co by nám v QGIS trvalo klikat hodiny můžeme v R naprogramovat během minut

# načteme všechny cesty k shp vrstvám v daném adresáři
parky.files<-list.files(path=cesta,pattern="*.shp$",full.names = T)
parky.files<-parky.files[-1]

## když chceme všechny soubor do jednoho objektu
# načte první shpfile
parky.init<-st_read(parky.files[1],quiet = TRUE)

# připojí všechny další shp file do jendoho objektu  
for (i in 2:length(parky.files)) {
  parky.next<-st_read(parky.files[i],quiet = TRUE)
  parky.init<-rbind(parky.init,parky.next)
}  

head(parky.init)
## Simple feature collection with 6 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 13.6537 ymin: 48.81219 xmax: 18.74347 ymax: 50.67169
## Geodetic CRS:  WGS 84
##   OBJECTID                       geometry
## 1     3900 POLYGON ((18.53572 49.66022...
## 2     3894 POLYGON ((17.88455 49.15734...
## 3     3873 POLYGON ((14.82518 49.66756...
## 4     3877 POLYGON ((14.19356 49.00166...
## 5     3904 POLYGON ((13.93738 49.81539...
## 6     3888 POLYGON ((16.10584 50.662, ...
# připojit metadat GIS join
meta.path<-paste0(cesta,"metadata_parky.csv")
metadata<-read.table(meta.path,sep=";",header=T,stringsAsFactors = T)
parky<-merge(parky.init,metadata)

head(parky)
## Simple feature collection with 6 features and 10 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 13.61426 ymin: 48.80132 xmax: 15.00087 ymax: 50.67842
## Geodetic CRS:  WGS 84
##   OBJECTID KOD  KAT                    NAZEV      ROZL  ZMENA_G  ZMENA_T
## 1     3873  21 CHKO                   Blaník  4029.195 20030101 20170901
## 2     3874  22 CHKO               Český kras 13225.701 20120824 20170901
## 3     3875  23 CHKO Kokořínsko - Máchův kraj 41037.143 20150223 20170901
## 4     3876  24 CHKO             Křivoklátsko 62497.146 20131107 20170901
## 5     3877  31 CHKO              Blanský les 21968.998 20091120 20170901
## 6     3878  32 CHKO                Třeboňsko 68744.620 20091218 20170901
##   PREKRYV SHAPEAREA  SHAPELEN                       geometry
## 1       0  40291954  31958.58 POLYGON ((14.82518 49.66756...
## 2       0 132257006  83895.23 POLYGON ((14.3244 50.00705,...
## 3       0 410371433 198507.51 MULTIPOLYGON (((14.4139 50....
## 4       0 624971460 138706.87 POLYGON ((13.81561 50.14976...
## 5       0 219689976  79790.35 POLYGON ((14.19356 49.00166...
## 6       0 687446196 143605.35 POLYGON ((14.71406 49.18412...
# kontrolní obrázek
plot(st_geometry(hrcr))
plot(st_geometry(parky),add=T,col="green")

# interaktivně
distCol<- colorFactor(distinctColorPalette(32), parky.init$OBJECTID)
leaflet(parky) %>% 
  addTiles() %>%  
  addPolygons(color = ~distCol(OBJECTID),
              stroke = FALSE, fillOpacity = 0.8,
              label = paste0(parky$KAT,"_",parky$NAZEV))

Rastrové analýzy

Velice dobrá knihovna pro rastrové analýzy v R je whitebox, což je implementace jinak v Python/Rust napsaných geo algoritmů.

## instalace balíčku whitebox pro rastrové analýzy
# install.packages("whitebox", repos="http://R-Forge.R-project.org") # instalace
# whitebox::wbt_init() # inicializace
library(whitebox)

# pozor na diakritiku a mezery v cestě. Pro WBT nesmí být
# cesta k rastrovým datům
dem <-("d:/Git/gisproba/data/DEM_Jested_100m.tif")

## stínovaný reliéf
# cesta k budoucímu výsledku
output<- ("c:/Users/Public/Documents/Hillshade_100m.tif")

# spustí algoitmus pro výpočet stinovaného reliéfu 
wbt_hillshade(dem, output, azimuth = 315, altitude = 30, 
                    zfactor = 1,verbose_mode = FALSE)
## [1] "hillshade - Elapsed Time (excluding I/O): 0.1s"
# načtu výsledek jako rastr
hill<-raster(output)
# načtu původní dem jako rastr
elev<-raster(dem)
# vytisknu oba obrázky pro kontrolu
plot(hill, col=grey(0:100/100), legend=FALSE,main="stínovaný reliéf")
plot(elev, col=rainbow(25, alpha=0.35), add=TRUE)

## Sklon svahu
output<- ("c:/Users/Public/Documents/slope_100m.tif")
wbt_slope(dem,output,zfactor = 1, verbose_mode = FALSE)
## [1] "slope - Elapsed Time (excluding I/O): 0.1s"
slp<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="Sklon svahu")
plot(slp, col=heat.colors(25, alpha=0.2), add=TRUE)

## The terrain ruggedness index (TRI)
output<- ("c:/Users/Public/Documents/TRI_100m.tif")
wbt_ruggedness_index(dem,output,zfactor = 1, verbose_mode = FALSE)
## [1] "ruggedness_index - Elapsed Time (excluding I/O): 0.1s"
tri<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="terrain ruggedness index")
plot(tri, col=cm.colors(25, alpha=0.3), add=TRUE)